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We investigate the interacting Dirac fermions on honeycomb lattice by cluster dynamical mean- 
field theory (CDMFT) combined with continuous time quantum Monte Carlo simulation (CTQMC). 
A novel scenario for the semimetal-Mott insulator transition of the interacting Dirac fermions is 
found beyond the previous DMFT studies. We demonstrate that the non-local spatial correlations 
play a vital role in the Mott transition on the honeycomb lattice. We also elaborate the experimental 
protocol to observe this phase transition by the ultracold atoms on optical honeycomb lattice. 
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The honeycomb lattice systems have been intensively 
studied recently, where a number of exotic phenomena 
are found in both experiment and theory, such as the 
single-layer graphene [lHSf, the topological Mott insula- 
tor |4[ and the recently found quantum spin liquid |5j. 
The electrons on honeycomb lattice can be described by 
a "relativistic" massless Dirac- fermion model [6[, which 
is characterized by the linear low-energy dispersion rela- 
tion E(k) = ±^i?|k|. For half-filling case, the Fermi sur- 
face of honeycomb lattice is reduced to isolated points 
at the corners of Brillouin zone, hence the density of 
states vanishes at the Fermi level, leading to the so called 
semimetal. When the repulsive interaction between the 
fermions is turned on, one can expect this semimetal to 
be driven into a gapped insulator by the sufficient strong 
interaction, in other words, a semimetal-Mott insulator 
transition happens. The Dirac fermion on honeycomb 
lattice is massless, hence the prevalent Fermi liquid the- 
ory is not proper here. Moreover, the semimetal has a 
small spectral weight near the Fermi surface, thus the 
low-energy many-body scattering effects would be quite 
different from those of the usual metallic systems. As a 
consequence of these anomalous natures, a new scenario 
is desirable for the description of Mott transition of the 
interacting Dirac fermions. There have been various the- 
oretical studies of phase transitions on the honeycomb 
lattice, however up to date no comprehensive conclusion 
has been achieved [7H9j. 

Due to the efficient description of the quantum fluctu- 
ations, the DMFT [10l4l3| and its cluster extensions [14J — 
llq have made substantial progress in the field of Mott 
transition in the past decades. The recent DMFT stud- 
ies on honeycomb lattice find a first-order Mott transition 
and the renormalization of the Fermi velocity vp induced 
by the interaction 0, [Hj|. Nevertheless, the DMFT ig- 
nores the non-local correlations. It is efficient only in 
physical systems with high dimensions or large coordi- 
nation numbers [To|. The CDMFT method [14], as a 
cluster extension of DMFT, effectually incorporates the 
spatial correlations by mapping the lattice problem into 



a self-consistently embedded cluster rather than a single 
site in DMFT. For low-dimensional systems, the quan- 
tum fluctuations are much stronger than in the higher 
dimensions, hence the CDMFT method is more precise 
than the DMFT here. To our knowledge, the DMFT 
calculations on square lattice basically agree with the 
CDMFT. The local approximation of DMFT captures 
well the main characteristics of the first-order Mott tran- 
sition [19j. However, in the ID chain lattice, the igno- 
rance of non-local correlations within DMFT leads to re- 
sults with severe errors [20|. The coordination number 
of honeycomb lattice is three which is between those of 
2D square lattice and ID chain lattice. Therefore, it 
is interesting to see whether the non-local correlations 
can significantly affect the behaviors of the electrons on 
honeycomb lattice and whether the CDMFT can provide 
novel results beyond the DMFT. 

In this Letter, we combine the CDMFT method with 



CTQMC simulation [2l|, which is employed as impurity 
solver, to investigate the interacting Dirac fermions. A 
second order Mott transition is suggested. We also find 
a number of novel features of the Dirac fermions, includ- 
ing the quasi-particles with anomalous long lifetime and 
the invariable density of low-energy excitations which is 
highly relevant to the non-local correlations. Our results 
are distinct from those of the previous DMFT studies 
[171 . Il8j , which indicate the non-local correlations are in- 
deed crucial in the honeycomb lattice. 

We consider the standard Hubbard model on the hon- 
eycomb lattice defined by the Hamiltonian, 
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where cr io . and q ct are the creation and the annihila- 
tion operator of fermions with site index i and spin 
index a, n^ = cl a Ci a is the density operator. U is 
the on-site repulsion, and \i is the chemical potential. 
The nearest neighbor hopping amplitude t (t > 0) de- 
termines the bare dispersion of the honeycomb lattice 

E(k) = ±tJ3 + 2cos(k y a) + Acos(y/3k x a/2)cos(k y a/2), 
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FIG. 1: (Color online), a) Illustration of the honeycomb lat- 
tice. The dashed line sketches the 6-site cluster scheme we 
explore in this work, b) The first Brillouin zone of the hon- 
eycomb lattice. The linear low-energy dispersion relation dis- 
plays conical shapes near the Fermi level. 




where a is the distance between the two adjacent sites, 
the plus and minus sign respectively denotes the up- 
per (tt) and lower (n*) band. At low energies, the dis- 
persion relation can be approximated as a linear one 
E(k) = ±u_p|k|, which can be described by the Dirac 
cones as shown in Fig. lb. We take t as the energy unit 
throught out this paper. 

The CDMFT method has been successfully applied to 
the studies on Mott transition of the Hubbard model for 
various systems [H-Q- We use CDMFT to map the 
original honeycomb lattice onto a 6-site effective clus- 
ter embedded in a self-consistent medium. The effec- 
tive cluster model is obtained via an iterative procedure 
which can be started with a initial guess of the cluster 
self-energy T,(zuj). The effective medium represented by 
the Weiss function g(iu) is determined by the cluster self- 
energy T,(iuj) via the coarse-grained Dyson equation, 
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where t(K) is the Fourier-transformed 6x6 hopping 
matrix with wavevector K in the cluster reduced Bril- 
louin zone of the superlattice, and \i is the chemical 
potential. Since the effective cluster model has been 
defined by a(iuj), we can employ the numerical meth- 
ods [2l|, I25U27J as impurity solver to calculate the clus- 
ter Green's function G(iuS). By using Dyson equation 
T^(iuo)=g~ 1 (iuo) — G -1 (zcj), we then recalculate the clus- 
ter self-energy Y,(ioj) to close the self-consistent iterative 
loop. This CDMFT loop is repeatedly iterated until the 
numerical convergence has been achieved. In the present 
work, we use the numerically exact CTQMC simulation 
as impurity solver and take about 7 x 10 7 QMC sweeps 
for each CDMFT loop. 

We first study the density of states (DOS) derived from 
the imaginary time Green's function G(r) by using max- 
imum entropy method [28|. The DOS for several values 
of on-site interaction U at temperature T/t = 0.05 are 
presented in Fig. 2. When the interaction is increased, 
the two quasi-particle peaks above and blow Fermi level 



FIG. 2: (Color online). The density of states for different 
interaction strength U at temperature T/t = 0.05. a) Within 
6-site CDMFT, the DOS near Fermi level is invariable until 
the Mott transition happens, which denotes a constant Fermi 
velocity vf in the semimetal regime, b) The 8-site cluster 
results. Insert: geometry of the 8-site cluster. Left panel: 
the low-energy DOS of central site A is invariable with the 
increasing of U, which agrees with the 6-site CDMFT result. 
Right panel: the low-energy DOS of boundary site B increases 
as U increases, resembling the DMFT result. Note that the 
Mott transition critical value U c of the 8-site cluster is much 
larger than that of the 6-site cluster and close to the DMFT 
result. This is due to the boundary effects. 



shift to the Fermi surface. The spectral weight is contin- 
uously transferred to the higher energy states, and even- 
tually a gap opens at the Mott transition critical value 
U c /t ~ 3.7. However, it is surprisedly noted that the 
DOS near Fermi level is independent of the interaction 
strength U until the Mott transition happens. The quasi- 
particles in semimetal regime are always analogous to the 
non-interacting Dirac gas, indicating no renormalization 
of the Fermi velocity. This result is completely different 
from the previous DMFT studies [17|, LL8| , where the re- 
sults suggest the DOS near Fermi level increases as U 
increases in the semimetal region, hence the Fermi ve- 
locity vf is renormalized by the on-site interaction. We 
argue that the non-local correlations which are absent in 
DMFT, but included effectively by CDMFT significantly 
reduce the low-energy spectral weight. 

To clarify the spectral reduction effect of the non-local 
correlations, we study an 8-site cluster which has two 
boundary sites with only one neighbor site in the clus- 
ter. The geometry of the 8-site cluster is shown in in- 
serts of Fig. 2b. It is apparent that the reduction of 
the coordination number subdues the non-local correla- 
tions for boundary sites [29j, therefore the characters of 
the boundary sites may be distinguished from those of 
the central sites. As shown in Fig. 2b, the central site 
shows invariable low-energy DOS coinciding with the 6- 
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FIG. 3: (Color online). Double occupancy D occ as a func- 
tion of interaction U at different temperatures. The double 
occupancy is insensitive to temperature for both weak- and 
strong-coupling regime, hence the three curves for tempera- 
ture T/t — 0.02, 0.05, 0.10 strictly superpose each other. The 
double occupancies decrease smoothly as U increases, indi- 
cating second-order Mott transitions. The transition points 
for different temperatures are marked by the arrows. 



FIG. 4: (Color online) The spectral functions A(k, uj) for a) K 
point. We can see the pseudogap at Fermi level developing as 
the interaction strength U increases, b) M point. The peaks 
of A(k, uj) are very sharp even in the strong coupling regime, 
which indicates well-defined quasi-particles. c) Y point. Far 
away from the Fermi level, the quasi-particle peaks are rapidly 
broadened by the interaction. 



site CDMFT study. However, the result of boundary 
site suggests low-energy DOS increasing as U increases, 
which is just a characteristic observed by the DMFT. 
The central and boundary sites are studied in one clus- 
ter, this disparity can only be attributed to the non-local 
correlations. We therefore conclude that the DMFT in- 
correctly suggests an increasing low-energy DOS due to 
the ignorance of the non-local correlations. 

We further investigate the double occupancy defined 
by the first derivative of the free energy F, D occ = ^ — 
(n^Hii) as a function of U for different temperatures. 
As shown in Fig. 3, we find that when the tempera- 
ture is lower than T c /t ~ 0.1, the D occ is independent 
of T in both semimetallic and insulating phase. Unlike 
the DMFT result [lij or CDMFT studies on triangular 
lattice [22[ and kagome lattice |24[ , where the double oc- 
cupancies D occ are insensitive to the temperature in only 
insulating phase. It is also remarkable that when U in- 
creases, the D occ curves decreases smoothly on the entire 
domain of U which indicates second-order phase transi- 
tions. There is no hysteresis phenomenon observed in our 
calculations, which also signals continuous phase transi- 
tions. This feature is expected, considering the fact that 
regardless of the on-site interaction strength, the upper 
and lower Hubbard band of honeycomb lattice always 
touch each other in only Dirac points in the semimetal- 
lic phase. Therefore the order parameter D occ decreases 
continuously even when the gap between the two Hub- 
bard bands opens. 

Fig. 4 shows the spectral functions A(k,(j) at K point 
k = (0,47r/3a) , M point k = (2 v / 37r/3a, 0) and T point 
k = (0,0). For the K point (Dirac point), as the in- 
teraction U is increased, U/t = 3.0 for example, the 
quasi-particle peak at Fermi surface (uj = 0) is strongly 
suppressed by the antiferromagnetic correlations, and a 
pseudogap develops. When U reaches the critical value 



U c /t ~ 3.7 for the phase transition, the pseudogap even- 
tually engulfs the whole Fermi surface and the system 
turns into Mott insulator. The observation of pseudo- 
gap on the Dirac point in the transition regime indicates 
the antiferromagnetic correlations are the origin of the 
Mott phase transition. Note that while the quasi-particle 
peaks at T point are rapidly broadened by the on-site in- 
teraction, the spectral functions A(k,uu) at the K and 
M point have sharp peaks even in the strongly coupling 
regime, indicating a long quasi-particle lifetime. This is 
related to the fact that there are few electrons near the 
Dirac points, thus the many-particle scattering effects are 
not prominent for the K and M points. 

Although the local density of states p(u) = ^2 k A(k, uj) 
is independent of the interaction for small energy uj, the 
states are redistributing in the momentum space when 
the interaction increases. As shown in Fig. 5a, the dis- 
tribution of quasi-particles with energy uj/t = 0.4 near a 
Dirac point (2-s/3tt/3, 27r/3) becames isotropic in differ- 
ent directions when U increases and finally the spectra 
averagely distribute in the whole ring. The Fermi surface 
denoted by A(k, uj/t = 0) is also presented. The first Bril- 
louin zone of honeycomb lattice has three K and three 
K points located at the corners of the hexagon, as shown 
in Fig. lb. Due to the particle-hole symmetry, the six 
Dirac points are actually all equivalent in our study. As a 
result, the Mott phase transition happens isochronously 
on the six Dirac points. Fig. 5b also shows that the lo- 
cations of Dirac points are steady in momentum space 
with the increasing of interaction U. 

Finally, we present the phase diagram of the interact- 
ing Dirac fermions on honeycomb lattice. The diagram 
is separated into semimetallic and insulating phases by 
the second-order Mott transition line. The critical value 
U c decreases as temperature T decreases, and a finite 
U c /t ~ 3.3 for zero temperature phase transition is sug- 
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FIG. 5: (Color online) The distribution of A(k,cj) in mo- 
mentum space for different interactions, a) The distribution 
of A(k,uj/t = 0.4) near a Dirac point (2^/3^/3,2^/3) be- 
comes isotropic in different directions when U increases, b) 
A(k,cj/t = 0) depicts the Fermi surface of honeycomb lattice 
reduced to six Fermi points. 



gested, which is quite agree with the large scale quantum 
Monte carlo result U c /t ~ 3.5 |5[. Comparing to the 
DMFT result U c /t ~ 10 Hi|, the Mott transition criti- 
cal value U c within CDMFT is much smaller, and more 
reasonable according to the previous studies [7H9| . 

On the experiment side, due to the high controllability 
and clearness, the interacting ultracold atoms on the op- 
tical honeycomb lattice ma y b e an ideal tool for the ver- 
ification of our results [30-32]. The optical honeycomb 
lattice can be produced by three intersected standing- 
wave lasers. Each of the standing wave is formed by two 
detuned laser beams intersecting at an angle of 70.4°. 
The ultracold fermion atoms like the 40 if is then loaded 
in the optical lattice [33|. The on-site interaction be- 
tween the fermion atoms can be adjusted from zero to 
strong coupling limits by the Feshbach resonance. Us- 
ing absorption imagining technology, the crucial physical 
quantities such as the double occupancy and the Fermi 
surface can be obtained to analyze the system [34J. 

In summary, we study the Mott transition of the inter- 
acting Dirac fermions on honeycomb lattice by CDMFT 
method. It is found that the low-energy density of 
states is independent of the interaction strength in the 
semimetal regime, hence the Fermi velocity always equals 
to the one of non-interacting Dirac sea. As the interac- 
tion increases, the system enters Mott insulator via a 
second-order phase transition accompanied by a pseudo- 
gap on the Fermi surface. This suggests the Mott tran- 
sition is dominated by the Slater mechanism. It is also 
shown that the non-local correlations significantly reduce 
the low-energy density of states of the honeycomb lat- 
tice. We argue that because of lacking non-local correla- 
tions, the DMFT method is inadequate in the study of 
strongly correlated fermions on honeycomb lattice. Our 
study may contribute to the comprehension of the in- 
teraction driven Mott transition and the recently found 
quantum spin liquid on honeycomb lattice [5J. 




FIG. 6: (Color online) Phase diagram of the interacting 
Dirac fermions on honeycomb lattice. The line separating the 
semimetal and Mott insulator marks a sencond-order phase 
transition. 
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